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A layered system of charges with logarithmic interaction parallel to the layers and random dipoles 
in each layer is studied via a variational method and an energy rationale. These methods reproduce 
the known phase diagram for a single layer where charges unbind by increasing either temperature 
or disorder, as well as a freezing first order transition within the ordered phase. Increasing interlayer 
coupling leads to successive transitions in which charge rods correlated in N > 1 neighboring layers 
^s^j , are unbounded by weaker disorder. Increasing disorder leads to transitions between different N 

phases. The method is applied to flux lattices in layered superconductors in the limit of vanishing 
Josephson coupling. The unbinding charges are point defects in the flux lattice, i.e. vacancies or in- 
terstitials. We show that short range disorder generates random dipoles for these defects. We predict 
and accurately locate a disorder-induced defect-unbinding transition with loss of superconducting 
^3 ' order, upon increase of disorder. While N — 1 charges dominate for most system parameters, we 

propose that in multi-layer superconductors defect rods can be realized. 



PACS numbers: 64.60.Cn,74.25.Qt,74.78.Fk 



I. INTRODUCTION 



There is considerable current interest in topological phase transitions induced by quenched disorder, a problem 
relevant for numerous physical systems. Such transitions are likely to shape the phase diagram of type II supercon- 
ductors. It was proposed^ that the flux lattice (FL) remains topologically ordered in a Bragg glass (BrG) phase at 
low field, and becomes unstable to the proliferation of dislocations above some threshold disorder (or field). The 
increased effect of disorder may lead to increased critical current, this providing one scenario for the ubiquitous and 
controversial "second peak"^ line in the phase diagram. Another scenario was proposed recently^ and is based on 
a disorder-induced decoupling transition (DT) associated with the loss of superconducting order, responsible for a 
sharp drop in the FL tilt modulus. An important question then is whether this DT occurs before the BrG instability 
(i.e within the BrG phase) or whether both occur simultaneously. 

Theoretically, two types of phase transitions were shown to be specific for pure layered superconductors. The first 
is decouplingSi^ at which the Josephson coupling as well as the critical current between layers vanishes. The second 
' is the proliferation of point "pancake" vortices, vacancies and interstitials (VI) in the FL above a temperature T^ef 
which, above some field, is distinct from melting, as shown in the absence of Josephson coupling 5 . It is believed that 
this pure system topological transition merges with the decoupling transition^ as the bare Josephson coupling is 
increased, being two anisotropic limits of the same transition 9 . This transition induces a loss of superconducting order 
(parallel to the layers by VI and perpendicular to them by the layer decoupling) while the positional correlations of 
the pure flux lattice is maintained^. This transition has also been studied in both limits in presence of point impurity 
q , disordes&ii as well as columnar disorder—. In particular, we have recently demonstrated^ the existence of disorder- 
induced, VI unbinding transition with loss of superconductivity in 3-dimensional (3D) layered superconductors, which 
. , would be particularly relevant to many layered superconductors and multilayer system a 2 ! 13 . 

Topological phase transitions in two dimensional systems are conveniently studied using mapping onto Coulomb 
gases of charges interacting via a long range logarithmic potential. Studying general three dimensional systems, even 
for pure systems, is considerably more difficult. The limit of layered superconductors with magnetic coupling only, 
provides one rare example where the problem can be studied analytically in 3D in a controlled way. Indeed in this 
limit the problem amounts to coupled layers with 2D Coulomb interactions. In the presence of quenched disorder, the 
problem becomes quite subtle already in 2D because charges can freeze into inhomogeneous configurations. Progress 
was made recently and it was showniiii 5 ^^^ that quenched random dipoles lead to a phase transition, via proliferation 
of defects at a finite threshold value of disorder, even at temperature T = 0. New analytical methods, based on RG 
for the charge fugacity probability distribution, and mapping onto a solvable model of directed polymer on the Cayley 
tree were developped in 2Di£ii£. In a short account of the present workii we have extended some of these techniques to 
study the 3D system in presence of disorder. Although a complete RG study along the lines ofi£ is possible in principle, 
we have used simpler, and we believe largely equivalent, methods. The first is an energy rationale which generalizes 
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the Cayley tree mapping. Second, we have introduced^ a novel Gaussian variational method which incorporates the 
effect of the broad fugacity distribution, a feature previously revealed by the R G 15 i 17 . This method was also applied 
to the single and two layer case in a related work on a random Dirac model relevant to quantum Hall systems^. 

The aim of this paper is to present details of our previous Letter— focusing on two themes. First we consider a 
general Disordered Coupled Coulomb Gas (DCCG) model system defined by integer ±1 charges on M layers in which 
the interaction 

energy between two charges on layers n and n' is 2J„_„/ lnr with r the charge separation parallel to the layers; in 
addition the charges couple to quenched random dipoles. A general study of this system is performed both via an 
energy rationale and by a variational method, with consistent results. These methods are explained in detail. Second, 
we apply this study to various physical situations, mainly to layered superconductors in an external field. We justify, 
stating clearly the assumptions, that VI in the vortex lattice of layered superconductors with no Joscphson coupling 
and in the presence of pinning disorder can be described by the DCCG model with quenched random dipoles. 

In section II we present the DCCG model and its mapping to a sine-Gordon type problem. In section III we 
develop a, T — energy rationale by an approximate mapping to Cayley tree problem. For the one layer case we 
find the well known critical disorder value of a cr — 1/8 for the onset of VI. For the many layer case we find that 
as the anisotropy r\ = —Ji/Jq increases a cascade of phase transitions appear at which the number of correlated 
charges on N neighboring layers increases. These "rod" phases appear at an decreasing critical disorder value until 
at rj — > 1/2 wc find N — > oo and a cr — ► 0. In section IV we develop an efficient variational method which is tested on 
the one layer system, allowing for fugacity distributions, knownii to be important in 2D since disorder becomes broad 
at low temperature. We reproduce the phase boundary in disorder-temperature plane separating an ordered phase 
(bound charges) and a disordered (unbound charges, i.e. finite VI density); the critical disorder parameter at T = 
is <r cr — 1/8 is recovered. We also find a first order line within the ordered phase (seen in the dynamics studjiii) 
which becomes a crossover line in the disordered phase. In section V we extend our method to the 2-layer system and 
find for the anisotropy 77 a critical value T] c = 1 — 1/ \/2 above which the single layer type transition is preempted by 
a transition induced by bound states of two vortices on the two layers with a cr < 1/8, in agreement with the energy 
rationale of section III. However, in a limited range of 1 - l/\/2 < r/ < 1/3 we find coexistence with a two gap state, 
which is not captured by the energy rational in its simplest form, but does not change the value of a cr . Of course, 
all of these above results truly involve renormalized values of coupling J„ en and disorder a ren . Although we have 
not attempted a full RG study, one main additional effect of RG is simply to substitute bare by renormalized values 
accounting for screening effects, which on the basis of the two layer case can be assumed to be small for our present 
purpose (i.e. identifying transition lines at low temperature). 

In section VI we develop the effective theory of layered superconductors with magnetic coupling between layers, 
but without Josephson coupling. We show that point disorder for the FL leads to quenched dipoles for the VI, hence 
the DCCG problem. For typical layered superconductors we predict the one layer type transition with an effective 
disorder parameter. However, by increasing the separation between layers, as in multilayer systems2*A£ to exceed the 
lattice parameter of the FL, one may realize the new N > 1 rod phases. 



II. MODEL FOR DISORDERED LAYERED COULOMB GAS 



In this Section we define the model for M coupled layers of disordered Coulomb gases and also in terms of an 
equivalent sine-Gordon model. Consider n(r, I) integer charges on the l-th layer at position r within the layer. The 
two-dimensional (2D) position r is defined on a lattice of spacing £, which for the superconducting system is the 
coherence length. We study the Hamiltonian: 

r^r' 1,1' r,l r,l 

where E c is the core energy, accounting for short scale energies r < £. Charges on the same or different layers interact 
with a 2D Coulomb interaction, with 

Irl 2ir 
G(r)«| rHoo lnU G,»h)^ ( 2 ) 

with G(r) = / G q (l - e 4q r ) and J — J t^jj (on a square lattice G^ 1 = i[2 - cos(g x £) - cos(g y £)]). Neutrality 
is assumed in each layer. The disorder potential Vi(r) can be considered as due to random dipoles. A dipole has 



3 



a potential ~ 1/r or ~ 1/q in Fourier space; hence the disorder potential on the Z-th layer V;(r) has long range 
correlations: 

V5(q)VWq') «9-0 2A i ^(2^) 2 ^ 2 )(q + q') (3) 
(fHWW^W)) «|r-^Hoo 4A„, In Jf^I (4) 

where A;;/ > 0. This logarithmically correlated disorder is the one which exhibits a phase transition - other types of 
disorder with either weaker or stronger correlations result in either ordered or disordered phases, respectively hence 
no phase transition as function of the disorder strength. One simpler case, which we will study in details, is the case 
of uncorrelated disorder from plane to plane, namely A/;/ = (jJq5u>. In that case one has 

|r-r'' 



M(r) - V^t')} 2 «| r -r< Hco 4a J 2 In (5) 

It is clear that the model on a square lattice defined by its partition sum Z\ at t — X){ n ( r I)} e ~ /3H can a l so be seen as 
a neutral 2D Coulomb gas model for Af-component vector charges. A given configuration of charges is thus defined 
by a set of vector charges {n(r, l)}i=i t ..M on a 2D lattice. 
We define the Fourier transform: 

»(q,fc) =d^^n(r,0e jq - r+lfcc " (6) 

l r 

where d is the spacing between layers and in a continuum limit — > / dz and £ 2 ^2 ~ * J d 2 r. The inverse formula 
for the charge density (per unit area) is 

n(r,l)/e= I I n(q,k)e-^ r - ikdl (7) 

J k J q 

with / q = I ^£f*> k = 2nm/Md with m integer, -[Af/2] + 1 < m < [M/2], and J fc = ^ £ fc * at large M, 

We perform disorder averages via the replica method, i.e. from the replicated partition function Z m in the limit 
ttl — ► 0, disorder averaged correlations and free energy are obtained. For integer m we have 

{n«(r,l)} 

with /? = 1/T, which on a lattice is exactly a Mm-component 2D vector Coulomb gas with integer charges at each 
site r with integer entries n a (r, Z) at each a = 1, ..m, Z = 1, ..M. The replicated Hamiltonian is2i 

= _ £ ^ 0)i , 6 na(rj /)G ( r _ r > b (r', Z') + (3E C £ n 2 (r, Z) (9) 

#I«,»'6 = - p 2 A w (10) 

summation over repeated indices is assumed unless otherwise specified. 

For system which is cyclic and (statistically) translationally invariant in the z direction, i.e: 

J w = J\i_ v \ A w = A|,_,/| (11) 

it is convenient to work with a Fourier space version which reads: 

^ M =^2/ [n a (q 7 k)(G ) ab (ci,k)nt( q7 k)+PE c J2n 2 a (rJ) (12) 
M JkJ 1 r,l,a 
4n 

(G ) afc (q, k) - -=-[ff(A)* 0l 6 - /3 2 A(fc)] (13) 



For later convenience we have defined g(k) = (3J(k), J(k) = dj^i Ji exp(i/crfZ), with J; = J k J(k) exp —ikdl. Similarly 
A(fc) = dj^i ^i,o exp(iZcc?Z), i.e. for disorder uncorrelated between layers A(k) = da J 2 . 
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We proceed to define an equivalent sine-Gordon system. We first rewrite the logarithmic interaction by using a 
scalar field Xa( r , I), 

Z^=( TT e " iXa(r ' " <,(r,0 ^ £;o 2]--.'.«"" (r,0 } x (14) 

{n a (r,Z)} r,Z,a 

The average is done with the weight exp[— i J k J q Xa(q, fc)(G ( 7 1 )afc(q, ^)Xb( c l; &)]> performing this Gaussian average 
one readily recovers Eq. JSJ. The inverse of Eq. (|13[) is derived by the inversion formula (A5 a b + -B) -1 = (1/A)8 a b — 
B/[A(A + Bm)), which for m — > yields 

(^) a6(q ,*) = g[^ Bl6 + ^]. (15) 

The product in Eq. (|14(l at each lattice point r can be written as a sum of all ±1, values of n a (r, I), i.e. a sum on 
all integer vector n = {n a .;}; a = 1, ..m, / = 1, ..M 

Z^=(W+ E F[„] e ^»,«— (r ' () ]) x (16) 

r {,^0} 

where the fugacity is Y[n] = exp[— f3E c J2 a i n a i\- At this point we make an approximation of small fugacities Y[n] 
(dilute limit) and write the above as an exponent 

^HlME Y[n]e^ n ^% x . (17) 

r {„^0} 

This approximation neglects harmonics of exp[n • x], i.e. it neglects vector charges with entries \n a ,i\ > 1. These 
harmonics are irrelevant near the actual phase transition—. Here and below we define n • \ — ^2 la n a ,iXa{ r > 0- The 
result Eq. (|17f) can now be identified as the partition sum for a sine-Gordon type Hamiltonian, 

PHsg = ^ [ /'xa(q,fc)(G - 1 ) Qb (q,fc)x^(q,fc)-^^r[n]exp l n.x(r). (18) 

Z JkJq r n 

where Xa(q ; k) = £, 2 dJ2 r E; Xa( r > i)e ifcdz+lqr . We note that the + sign for the off diagonal replica term in Eq. <(T5|) 
corresponds to imaginary gauge disorder in a related Dirac probleml&. The validity of the approximations leading to 
(I18|l are discussed ini£ in the context of a single layer. As also shown below, it is important, as done here, to retain 
replica charges with several non zero entries in order to describe the freezing transitions at low temperatures. 



III. ENERGY RATIONALE 



In this section we consider the Coulomb gas problem at T — and develop an energy rationale to determine the 
phase diagram of the coupled layer system. The problem amounts to find minimal energy configurations of charges in a 
logarithmically correlated random potential. To ascertain the XY ordered phases (bound defects) and the transitions 
out of it (defect unbinding), a first step is to study the dilute limit of a single charge (or dipole). Even then, the full 
analytical solution is difficult, but various approximations have been argued to give exact leading order results. For a 
single layer it was studied either us i n giii2i22i^S a "random energy model" (REM) approximation, or more accurately 
using a representation in terms of directed polymers on a Cayley tree (DPCT), introduced iniSiSS. The continuous 
version of the DPCT representation (branching process) was shown to emerged from the one loop Coulomb gas RG 
of the single layer problem, both for the single charge (or dipole) problem and for the many charges problem including 
screening effects. It is thus expected to be accurate. 

Schematically, one considers a tree with independent random potentials (Fig. 1 inset) Vi on each bond with variance 
vf = 2cr Jq. For definitcness we can discuss a tree of coordination e 2 , the choice being immaterial for our present 
considerations. After p generations one has ~ e 2p sites which are mapped onto a 2D layer: each point r corresponds 
to a unique path on the tree with v±, v p potentials and is assigned a potential V(r) = v± + ... + v p . Two points r, r', 
separated by |r — r'| ~ e p in Euclidean space, have a have a common ancestor at the previous p' w In |r — r'| generation 
Since all bonds previous to the common ancestor are identical [V(r) — V(r')] 2 = 2^^ =1 u 2 = 4crjQln(|r — r'Q, 
reproducing Eq. (JSJ) on each layer. Thus the growth of correlations on the tree and in Euclidean space is by construction 



5 



the same, and the single charge problem corresponds to a single directed polymer. Exact solution of DPCT 20 yields 
the best energy gained from disorder V m i n = min r V(r) as —\/8aJo In L for a volume L 2 , with only 0(1) fluctuations^, 
i.e — V8aJo per generation p = InL. It is argued that this is also the exact result for the Euclidean problem. For 
a dipole in a single layer, one consider two directed polymers on the same Cayley tree. Opposite sign charge see 
opposite disorder — Vi, the gain from disorder —V max behaving identically. The configurations of the two oppositely 
charged polymers can however being argued to be essentially independent (i.e. determining maximum and minimum 
of a log-correlated landscape can be performed independently). 

To generalize the Cayley tree argument we construct optimal energy charge configurations for M coupled layers as 
follows. Consider N neighboring layers with a +, — pair on each layer and no charges on the other M — N layers. 
We assume, for convenience, that Jo > and J;^o < so that equal charges on different layers attract. The DPCT 
representation now involves on a single tree N + polymers (each seeing different disorder) and N — polymers (each 
seeing opposite disorder — Vi to their + partner). A plausible configuration is that the + charges bind within a scale 
L e (0 < e < 1), so do the — charges, while the + to — charge separations define the scale L. Its tree representation 
(Fig. la) has 27V branches with elnL generations, i.e. an optimal energy of — 2/Vv / 8crJoem£. On the scale between 
L e and L the + charges act as a single charge with a potential X)z=i ^i( r ) (t ne N polymers share the same branch) 
of variance Na hence the optimal energy is — 2V8iVcrJo(l — e) InL. Note that the rod formation limits the disorder 
optimization leading to a disorder energy ~ y/N < N. The total energy gain from the disorder potential is thus 
estimated as: 

E dls w -2J %/8^[eiV+ (1 - e)y/N]]xiL. (19) 

It is clearly exact for both e = and e = 1, sufficient for our purpose. This result can also be obtained from 
the REM approximation, i.e. replacing the V(r) by L 2 variables uncorrelated in r, with the same on-site variance 
V 2 (r) ~ 2a Jl In L also yielding^ V mm VSa J In L. 

The competing interaction energy from the couplings Ji is for the H — pairs [2JqN + 4^^ Ji(N — I)] lni while 
for the ++ and pairs it is — 4^^ 1 Ji(N — l)e\nL. Hence the interaction energy is 

N 

E int = 2J N[l-2j2 m (l-l/N)(l-e)]lnL (20) 
i=i 

where rji — —Ji/Jq. The total energy Etot = Edis + -EWit is linear in e, hence the minimum is at either e = 1 or at 
e = 0. Since e = 1 implies that the + charges unbind, it is sufficient to consider e = with all N > 1, i.e. a rod is 
aligned with N correlated charges at distance O(l) and has energy 

N I /8a 

Etot = 2J N[1 - 2 W(l - jj) - J ^ ] InL . (21) 

i=i 

One can introduce more scales L e to describe the multi charges, however, as the energy is linear in e' the result is 
the same rod structure. 

Consider first the case with only nearest neighbor coupling 771 and only intralayer disorder correlation <jJq . Disorder 
induces the N vortex state (i.e. E to t vanishes) at the critical value 

^ ) = f[l-2r 7l (l-i)] 2 . (22) 

The system is thus fully stable to disorder only for a < a cr with: 

. N, „ ,„ 1 

run — 1 — 2r?i (1 

n 8 L ln N 



mm^[l-27 ?1 (l-^-)] 2 . (23) 



When a reaches a cr the first instability is to one of the N rod state, where iV depends on the value of the anisotropy 
771. If t]i < r/^ = 1 — l/\/2 then air is minimal at N = 1 and the first instability is similar to the one of a single 
layer with a cr = 1/8. For larger anisotropies ri[ N < i]i < r/[ N \ the first instability occurs at a cr = a c ^ towards 
a 7V_ r od state with 1/(1 - 2r/[ N) ) = 1 + ^N(N + 1) - N, thus with diverging N as r][ N) -> \ (Fig 1) (for 771 > \ 
Etot < even without disorder and the defects would form a lattice at T = 0). 

Upon increasing a beyond o'er a given rod phase N > 1 would eventually decompose into the N = 1 phase. 
In particular the energies of the N = 1,2 phases become equal at a c l' 2) = ?7 2 /[4(V2- l) 2 ] which equals 1/8 at 
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FIG. 1: Critical disorder values with only nearest neighbor coupling Ji vs. the anisotropy r\ — —Ji/Jq. Transitions between 
different N phases are marked with arrows. Inset: the Cayley tree representation (for iV = 3 neighboring layers) with + charges 
(at the tree endpoints) separated by L e along the layers, and separated by L from the N — 3 ~ charges. 



7]i = 1 — l/\/2- Hence at rji > 1 — l/\/2 the N=2 rods disintegrate into N=l charges at a > a^r' 2 \ The variational 

( 1 2) 

solution (section VB) shows that this secondary line is actually at a somewhat lower a cr ' ' (see Fig. 5 below). 
In the general case with all couplings J; the critical value is: 

aW = ^[l- 2 f>(l- I)] 2 . (24) 
i=i 

We consider in particular J; with range of Iq constrained by J; = 0, as relevant for the superconductor system 
(section VI). An illustrative example is rji — 771 exp(— (I — 1)/Iq), constrained as 771 = |(1 — exp(— 1/Zq)) (note that 

YI1L1 Vi = 5(1 — ex P( — N/lo))- One then has: 

{N} 1 (1 - exp(~iV/Z )) 2 
acr 8N (1 - exp(-l// ))2 1 ' 

For large Iq 1, each rji^o is small: function of AT starts by increasing and for N < Zo the lowest cr is at 

AT = 1. However, the combined strength of N « vortices being significant, it has a maximum and then decreases 
back to zero for N > Iq as 0-^' « Iq/8N. Hence air — > as A — > 00 and any small disorder seems to nucleate such 
vortices. This is because of the perfect screening of the zero mode J; = which implies that an infinite charge rod 
has a vanishing lnr interaction; hence a logarithmically correlated disorder is always dominant. 

In practice, the realization of these large N states depends, however, on the type of thermodynamic limit. Adding 
to Eq. if^Tl) the core energy 2E C N yields 



E' tot = 2J Q V~N(y 8a%P - V8a) ln£ + 2E C N (26) 
which becomes negative only beyond the scale 



L N nexp{E c VN/[Jo(V8^-^8-a£>)]}. (27) 

This Ljv is the typical distance between rod vortices. Hence even if a > air only for system size L > Ln the energy 
gain from disorder wins over core (and interaction) energy. Hence as a — > such states are only achievable in a 

thermodynamic limit where L/N diverges exponentially. Using air ~ Iq/SN, for A > Iq/8<j the lowest scale L in 
this range is achieved at A — Iq/2& and leads to a (system size) lower bound L m i n w exp[E c lo/4Joa] for observing 
large A states with a given a < |. For layered superconductors^ E c /J ^S> 1 and l ^> 1 and this large A instability 
occurs at unattainable scales, thus A = 1 dominates. One needs lo as 2 — 3 to realize the A > 1 states, attainable in 
multilayers (see discussion section VII). 

We finally generalize the energy argument for the e = configuration to the case of arbitrary correlations 7; = 
Aj/Aq. The disorder energy can be found from the variance of X^ijv^K 1 ") leading to the replacement a — > 
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cr(l + 2 ^ i=1 (A; i o/Ao)(l — l/N)) in Eq. (|21|l . A more compact form can be obtained by writing directly 

N N 
Etat = E Jl1 



1 = 1 l'=l 

Using that: 

N N 



N N 



\ 8 EE A "'i lnL ( 28 ) 
\ i=i i=i 



EE J «' = [ J(k)4> N (k) (29) 
i=i i'=i Jk 

N (k)=j:e^ = ^^§ (30) 
'—' sin (kd/2) 



w 

One has the criticality condition for a N rod: 

A(fc)<Mfc) = i( f J(k)<f> N (k)) 2 (31) 

'fc s Jk 

which in terms of a = Aq/Jq has the critical value 

(JV ) 1 Ufe ggTgTg ■W-*)) 

^ ~ 8 J.^rSTf A (fc )/A ( } 

For fixed anisotropics J(k)/Jo, A(fc)/A this relates the overall critical disorder strength to the rod length N . 

IV. VARIATIONAL METHOD - THE SINGLE LAYER 

We develop here a variational method which allows for fugacity distributions, an essential feature in the one-layer 
problem. The method is developed in this section for the one-layer system and it is shown that one recovers in a 
simple way all the important known features for this problem. Furthermore, new insight is gained for a critical line 
within the ordered (charge bound) phase, as well as a crossover line in the disordered (charge unbound) phase, at 
which the the functional dependence of the charge density changes. 

The single layer replicated Coulomb gas Hamiltonian is 

f3H {m) = \f ^(q)^[^- ff A' 2 K(q) + /JE c ^^(r). (33) 
1 J i q 

where n a (q) = J2 r n a( r ) e ' q r - Note that a > is here essential; the same 2D system with a < has been shown to 
have a different phase diagram2iS4. The equivalent sine-Gordon system is now 

f3HsG = \ I ^(°l)( G o 1 )a6(9)X6(q)-EE y [ n ] ex P in ^( r ) ( 34 ) 

•'I r n 

(Go 1 )a b (q) = ^[p ab + ^} (35) 

where Xa(q) = £ 2 5Z r Xa(r)e lq ' r and bare fugacities Y[n] = exp(— (3E C J2 a n a)- Here one has simply n • x( r ) = 
n aXa(r), with n a nonzero vector with entries n a — ±1, 0. 
The variational method represents the full Hamiltonian (|34ll by an optimal Gaussian one of the form 

mvar = \j G- b 1 (q) X a(ci)X b (~<l) (36) 

where G a b is to be determined by a variational principle. The variational free energy is F var = Fq+ < Hsg — 
Hvar >H % , ar with (3Fq = — IiiZq — — \Tr In G is found to read: 



PF va 
L 2 



-\ f Tr\nG{q)+ 1 - f T^Gq 1 (q)G(q)) - £~ 2 ^ Y[n]e ^ J* n ' G(9) ' n (37 ) 
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up to an unimportant constant, where the Tr is in replica indices. Taking the derivative SG one obtains the 
saddle point equation: 

where we have defined: 

G-J{q) = {G Q )-J(q)+cx ab (39) 

We recall first some technical relationship. In the following we represent relevant operators as averages which 
depend only on n + ,rt_, which are the number of + or — entries in n, respectively. The averages have the form 

A[n] = J dudv^{u,v)e u{n + +n ^+< n +- n -^ (40) 

where z± — e u±v can be interpreted as fugacities for the ± charges, hence $(it, v) is a fugacity distribution. A sum 
on all n =/= can be written in terms of the variables n + , n_ with a combinatorial factor for the number of n vectors 
with a given n + , n_ , 

li m J_ A[n] = lim — V — = r—^ / e («+^)n++(«-^)n-\ $ (4i) 

m—>o ui — ' m^o m ^— ' n+\n-'.(m — n + — n_ ! 

n^O 0<n++ri_<m ^ v ^ ' 

= lim < (1 + e u+v + e u - v ) m /m >$= (ln(l + e u+v + e u ~ v ))^ (42) 

?n—>0 

and the binomial expansion has been used and < ... >$ denotes an average with the weight $(u,w). Similarly one 
has: 



Hm n 1 E E = (E( n + + ™-)e" (n++ "- )+t 

m^O m. — ' ^ — ^ ^ — ' 



hv(n + -n_)\ 



= <9 U ln(l + + c— )). = ( 1 + ett+ T +eM _J , (43) 



and 



, - n _\2 e u(n + +n_)+»(n + -r._)\ 

m^O m ' — ' L — ' " " — ' 

n^O a,b n/0 



p u+v , u-v i 4„2ii 

= (5^ ln(l + e u+v + e u - v ))* = (- —-^)* (44) 

x v v n ( 1 -)- e ^, + ^, -|- e u ~^)2 

In our case we consider a replica symmetric parametrization a a b = <7 c 8 ao + (To so that G a b — J q G ao (q) has the form 
Gat = G c 5 ab - A, where 

f 1 A 2 
G c = — 2 = Kln (45) 

A=[ ao 2 + & = irVln-^--^V+^ (46) 

where A ~ ^S> Ka c is a cutoff on the q integration. Since J n • G(q) ■ n = G c n a ~ C(Sa n a) 2 we can now 
identify the weight function from the interaction term in Eq. (|37J) . 

Y[n]e~* X, n ' G( « } ' n = e -(iG c +^)(n + +„_)+ii(n + -n_) 2 

dv$(v)e u( "+ +,l - )+l,(,l +-™- ) (47) 
where here the weight function depends here only on v 
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and u = —f3E c — \G C . We recall that: 

y = e- 0Ec (49) 

is the bare fugacity of the charge, while the z± corresponds to the renormalized ones (they become random 
variables because of the quenched disorder in the system). The bare model can be generalized by introducing 
short-ranged randomness in the bare core energies (of width E )^H, resulting in the replica symmetric form 
Y[n] = exp[— j3E c J2 a n 1 ~ J2ab n a n b]- This corresponds to the change A — > A + fi 2 E 2 in the averages above. 

Since A is divergent at criticality a finite E can be ignored. 
The interaction term in Eq. I|37(l is therefore 



Y{n]e U, n ' G ^- n = (i n (i + e u +" +e u ~ v )) 9 (50) 



E 

To identify a c , oo we consider the variational equation (|38l) and note that Eq. I|44[) in the limit m — ► is tr c , while 
Eq. (|4*S|) is a c + <r a , hence 

p U+V _|_ U-V I A p 2u 

a c = r 2 ( ^ ! -^tt)* (51) 

- ^ + e u+v + e u-vy 

These equations, together with (|48(l and l|46() form the closed set of self-consistent equations that we want to solve. 
On general grounds one expects an ordered phase where the self energy a c vanish corresponding to zero charge density 
and zero renormalized fugacity (XY phase). The solution with a c > corresponds to a phase with finite density of 
charges (disordered phase), the typical correlation length (see Eq. I39fl being ~ cr^ 1 ^ 2 , the typicall distance between 
charges. We will thus perform the analysis near the critical line, where a c is small. We will first neglect the <jq term 
in H4tj|) and later show that it is indeed negligible in all regimes of interest. 

To analyze these equations we note that the v integration is dominated by large |u| and A which diverge at criticality, 
a c — * 0. The function displayed in l|51|) is maximal at v — — u with a width O(l), while the gaussian $(i>) is maximal 
at v — with a width 0(V~A). Consider then v > where the e u+v term dominates and is either very small (u + v < 0) 
or very large (u + v > 0) , hence 

J l + e u + v % /2^4 Jo V2^A J~ u V^A 

In the second term the saddle point at v = —A is outside of the integration range, hence it is dominated by the lower 
limit, i.e. it is of order exp(— u 2 /2A). The first term has a saddle point at v — A which is within the integration range 
if A < — u and then 

a c ~ e u+A ' 2 ~ yai K - K2 ^ 2 A < -u. (54) 

For A < —u the second term of l|53f) is indeed smaller, exp(— u 2 /2A) < exp(u) <§; exp(it + A/2). The range where 
a c is finite, i.e. the charge density is finite and behaves as a plasma is where the exponent in the solution is positive 
(both a c and y being small), 



<7 C - y2-K +K ^ k - K A a - 2 > 0, a < 2jK z (55) 

and the critical line where a c vanishes is K ~ K 2 a — 2 = (Fig. 2); the condition A < —u becomes a < 2/ K 2 (see 
below). This is the first, or high temperature regime. In that regime a standard small fugacity expansion works, the 
effects of the width of $(v) are unimportant, both at the transition and in the disordered phase. 

Considering now the second, or low temperature, regime A > — u. Then the first term of (|53|l is dominated by the 
upper limit, hence both terms of (j53|l yield 



e 



-u*/2A „ y i/*K„ a i/*, ff l A >-u. (56) 



Note that this corresponds to the distribution $(u) being very broad and then the maximum at v — — u dominates 
the result. For the finite charge density phase we have now 

<J C ~ y^ 1 a- > i cr > 2/K 2 (57) 



10 




o.oo 



0.25 



0.50 



T/J„ 



FIG. 2: Phase diagram for one layer in terms of a and T/Jo = 1/K variables. The full line is the defect transition given by 
K — K 2 a — 2 = 0ati<l/A"<i and by (7 = -| at 1/-BT < ^. The dashed line a = 2/K 2 is a first order transition within the 
ordered (low T) phase while a crossover line in the disordered phase. 



so that the critical line is a — | (Fig. 2); the condition A > —u becomes a > 2/K 2 (see below). 

The boundary between the regimes (|55|l and (|57|l is A = —u, which for cr c — > is a — 1/4K, i.e. a = |, K = 2 on 
the critical line. The form l|55|l is then valid at high temperatures K < 2 and a sum on single replica, single charge 
excitations is sufficient. In the low temperature regime (K > 2), where Eq. 1)57(1 is valid, the summation on all charges 
in all replicas n a = (0, ±1) is essential in obtaining the correct result. It corresponds to the physics of the freezing, or 
prefered localisation of the charges in deep minima of the random potential. 

It is instructive to evaluate the boundary between the regimes l|55[) and (|57|l for arbitrary small bare fugacity y <C 1 
also away from the critical line. The non-analytic behavior of the integral in Eq. 1)53(1 is related to the divergence of 
w, i.e. it exists in the ordered phase, while becomes a cross-over line in the disordered phase; this is further discussed 
below. Consider then a finite a c and define a c ~ y 7 . For y <C 1 the condition A = —u becomes a = (1/2K) + (1/jK 2 ). 
For A < —u we have from Eq. 7 = 2/(2 — K + crK 2 ), hence the boundary is a — 2/K 2 . Similarly, for A > —u, 
using u = (Kj/2 + 1) \ny yields 7 = (2 / ' K) j '(V8<t — 1) and again the boundary is at a = 2/K 2 . Hence there is a 
unique boundary between the two regimes [as included in the conditions for Eqs. and fBTjl ] which intersects the 
critical line at a = g, K = 2. The sharpness of this boundary, as mentioned above, depends on a c — > 0, hence in the 
disordered phase it depends on the smallness of y, i.e. it is a crossover line where the charge densities ~ a c change 
from Eq. (|55() to Eq. JSJJ, a crossover whose width shrinks with y. In the ordered phase a c = and formally the 
boundary is sharp, although the relevant observable, i.e. the density, vanishes. One may still observe this transition 
by a finite size effect where the q — > singularity is cutoff by the inverse area 1 / L 2 instead of a c , i.e. <r c ~ (1/ L) K ~ aK 
or ~ (l/L) 1 / 4a in the two regimes, respectively. This transition is termed as a freezing transition; it is related to the 
single directed polymer transition on a Cayley tree2£, to a dynamic transition^ and also to a phase transition in a 
random gauge Dirac systemic. 

Consider next ctq, Eq. (|52|l . The integral is again dominated by large |v|, hence 



CO 



- v 2 /2A dv 



(1 



V^kA 



2u+2v-v 2 12 A dv 



V2^A 



2 /2A_ 



dv 



^/2^^A 



(58) 



The second term is ~ exp(— u 2 /2A) while the first term has a saddle point at v = 2A which is inside the integration 
range if 2A < —u, and then 



a - y 2 e 2u+2A ~ y 2 (a c ) K -^ K a < 2/3K 2 . 



(59) 



2A < -u implies a < (1/4K) + 1/(2(3K 2 ) and since also A<-awe can use (3 = 2/(2 -K + a K 2 ), hence a < 2/3K 2 



Note that ctq ~ o". 



2-aK z 



<C cr c when a < 2/3K 2 , so that ao/<j c in Eq. (|46|l can be neglected. Consider next 2A > —u 



11 



where the integrals for oo are dominated by the end points v — —u. The range — u < 2A < —2u which corresponds 
to 2/3 A" 2 < cr < 2/K 2 yields 

cro ~ ( (7c )d+^ 2 /2) 2 /(2^ 2 ) 2/3K 2 <o< 2/K 2 (60) 

for which again <tq <C <J C while at a > 2/K 2 we have <7o ~ °"c- At cr = 2/3A 2 the functional form of Co changes, but 
since near this line Co <C a c there is no observable singularity. 

To conclude, comparison with RG studiesi^ii shows that the present variational method, which accounts for broad 
fugacity distributions, gives a remarkably accurate description of the transition and in particular of the freezing 
phenomena at low temperature in the single layer model. This is presumably because the screening effects (neglected 
in the variational approach) was shown, via higher order RG, to be very small at low temperature. In addition it 
provides a description of the disordered phase. The RG methods can be extended to many layers but following the 
joint distribution of the large set of random fugacities becomes rapidly difficult as M increases. We thus now turn to 
the extension of the variational method to several layers. 



V. VARIATIONAL METHOD - MANY LAYERS 



A. general case 

We study now the full many-layer system, Eq. (|18fl . We develop a variational method for M coupled layers which 
allows for fugacity distributions, an essential feature in the one-layer problem. It is explicitly worked out for two 
layers, describing the various rod transitions as found by the energy rationale in section III, as well as a narrow 
transition region. 

We note in particular the form of the interaction term expin • x( r ); the naive approach would be to restrict to 
charges n with a single non zero entry, leading to a uniform fugacity term -J/D rlla cos(x na (r)) and a diagonal 
fc-independent replica mass term. Instead we keep all composite charges n, which allow for variational solutions with 
off diagonal and fc-dependent replica mass terms. This corresponds respectively to fluctuations of fugacity and N > 1 
charge rods being generated and becoming relevant. 

We note first that a rod solution is readily obtained from Eq. i|12[) . i- e - we l°°k f° r N correlated charge on nearest 
layers so that 

N-l 

\n a (q,k)\ 2 = K(q)d^ e ife T = d 2 K(q)| 2 <Mfc) (61) 
;=o 

where </>jv(fc) was defined in l|3U|) . With this replacement Eq. 112fl has the form of a one-layer system Eq. (!.').' jfl with 
effective parameters 

Ke ff = [ g(k)<p N (k) (62) 

= (JJik^MW (63) 

The system than has the same phase diagram as for one layer (Fig. 2) with these effective parameters. In particular 
the T = transition is at <J e ff — \i hi agreement with Eq. (|31|l . 

We proceed with the variational scheme and define an optimal Gaussian Hamiltonian to approximate Eq. 118|) as: 

Hq = \JJ Xa(q,A)G„ 6 (q,fc)xJ(q,fc) (64) 

G- b \d,k)) = (Go)- b 1 (q,fc)-l-a c (fc) ( 5 afc + cr (fc) (65) 

i.e the self energy can now depend on k. 

The variational free energy is T va r = + (Ti-sa — ^o)o where (...) is an average with respect to Ho and JFo is its 
free energy f3Fo = — \Tr In G a fc(q, k). The Gaussian average has the form 

F[n]^(exp*n. X (r)) = cxp{-i / ^ | ]T m, a e Ml \ 2 G c (k) + ~ / J £ n^e lkdl \ 2 A{k)} (66) 
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where we recall that J k = J2k and I q G ab(q, k) = G c (k)S ab - A(k), with 
G c (k)= [ G c (q,k) = .g(fc)ln[A 2 /(4ng(k)a c (k))] 

Jg 

A(k) = [ G(q, k) = (3 2 A(k) ln[A 2 /(4^ 9 (fcK(fc))] - P 2 A(k) + g(k)a (k)/a c (k) . (67) 

Jq 

Extremization of F var yields the saddle point equation: 

{vch'Sab + Mi' = r 2 d~ 2 E n al n w Y[n]F[n] (68) 

n 

since the dependence is on I — Z', a corresponding Fourier transform yields 

<?c(k)6 ab + a (k) = r'd" 1 E E n al n bl ,e ikd «- l ')Y{n}F[n} (69) 

n l-l> 

We can now define s a (k) — d^ jl n^ a e lkdl . The A(k) term can be written as an average over gaussian distributions 
of fugacities: 

n «p (^i-a(*)i^(*)) - J n ^ «p e - a^Ki a ) ( 7 °) 

This form allows to decouple J2 S F[ s ] = (Z m )to with 

Z = E ex P(-^ / Gc(fc)| S (fc)| 2 + - f H[ Wfc **(fc)]) (72) 

{s„=0,±l} k k 

The variational equations for m — > become 

'•<*>- f^!s5|>" ! "«w = «- 2 "<i^fi 2 >« <™> 

We will not attempt to solve the general case but rather present a solution for M = 2. 

B. detailed solution for two layers 

We consider now two layers with uncorrelated and equal disorder on each layer. The partition sum depends now on 
the number of + and — charges on each layer, i.e. on the 8 numbers n a ,/3 where a, (3 = ±1,0, excluding noo- For the 
vectors ni, n 2 in replica space for each layer, their number for a given collection of n ai p is the combinatorial factor in 
the following sum 

E Y [m, n 2 ]F[ ni , n 2 ] = £ — — ^ exp[-/?E c ^Ki + n 2 2 ) - i-G c (0) 5> ol + ^2? 

ni ."2 not, (3 a a 

-^G c (tt)EKi - n a2 f + ^A 1 (J2n a i+n a2 ) 2 + ^A 2 (^n al - n a2 ) 2 ] (74) 

a ~ a * a 

where Ai = A(0)/2d, A 2 = A(n)/2d and the sum is restricted to ^2 a p n a ,p = m . We need then two fugacity 
distributions, 

exp[iA l( ^n al ±n a2 ) 2 ] = J ^^H^A (75) 
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where the upper (lower) signs corresponds to i = 1 or i = 2, respectively. The sum over n a ^ has the form of a 
"ninomial" expansion, i.e. a power of 9 terms, 

53 Y[n 1 ,n % ]F[n 1 ,n 2 ]={Z m ) u (76) 

ni,n 2 

where the average is on both <j\, u> 2 . In terms of u\ = —^f3E c — j^G c (0) and u 2 = — \fiE c — ~i^G c (%) we have 

^ _ J _)_ g"l + "2+^l+W2 _|_ g«l+tl2+^l— "2 I g«l+t'2 — Wi+U 2 _|_ pMl+M2-Wl-tJ2 

The equations for the (dimensionless) self mass terms a c \ = ^2^<7 C (0), a c2 — ^j^a^ir) and similarly for aoi are 

o C i =< 5 2 >^ (78) 
(dZ/d^) 2 

a- 0l =< 7^ >w ■ (79) 

These self masses correspond to length scales, i.e. a cl is the typical distance between (++) charge rods (i.e. a 

— 1/2 

+ in layer 2 is on top of a + in layer 1), while a c2 is the typical distance between (H — ) charge rods. In general 
Cc2 ~ [cci] Q so that a — corresponds to a c \ — with N=2 (H — ) rod defects, a = oo corresponds to er C 2 = with 
N=2 (++) rod defects, a — 1 corresponds to the two length scales being equal hence an N=l state, while other values 
of a imply the presence of two independent length scales. 

AJV = 2 rod solution is readily obtained by a c2 = so that u 2 -> oo and Z = 1 + e 4 " 1+2wi + e 4 «i- 2 "i. This is 
equivalent to the one layer system with 

K eff = 2(K ±K 1 ) 

*» = 2 { kI±K 1? (80) 

where the lower sign corresponds to the +, — rod solution a c \ = 0. 

Consider now a general solution of the form a c2 ~ [c c i] Q so that near criticality 



ui -> -i(JT + #i) ln(A 2 /a cl ) Ai -» iaif 2 ln(A 2 /a cl ) 



!X2^-Ja(^-^i)MA 2 /(7 c i) A 2 -> iaa^ 2 ln(A 2 /a cl ) (81) 

Near criticality the u> integrals are dominated by large values so that positive and negative integration ranges are 
equivalent; furthermore, the u>i,u> 2 > integral is dominated by exponents where lu\,lu 2 appear with positive sign, 

d dZ/duji d e «i+«2+^i+^2 + 2e 4tll+2a)1 

CTcl ~^a^ Z " 4 ^9^ 1 + e»i +«2+^i+^2 + g^i+^x + e 4„ 2 +2^ W=>o ( 82 ) 

We focus here on the low temperature behavior where — > oo and the integrals are dominated by the maxima of 
the above d/du>\. The fraction in Eq. (|82() has values 0, 1,2 as indicated in Fig. 3 with boundaries shown by the 
full lines, assuming for now u 2 > u\ (the solution for u 2 < u\ can be inferred by the symmetry of the phase diagram 
under Ko,Ki,a — > K\, — Kq, l/a). At the full lines in Fig. 3 d/du\ is maximal and dominate the integral at low 
temperatures since the Gaussian averaging factors are very flat. More precisely, we have separated the u)\ integral into 
ranges left and right of these lines and check in each range that it has no saddle point and is therefore dominated by 
u>i at the line position. Thus for u> 2 < —2u 2 the integral is dominated by lo\ = —uj 2 —u\ — u 2 leading to a contribution 

r- 2u 2 

while for lj 2 > — 2it2 the intgeral is dominated by u)\ — lo 2 + 3u 2 — u\ with the contribution 



4 2) ~ / duj2e -(^+3u 2 -u 1 ) 2 /2A 1 ~^ 2 /2A 2 (84) 
2u 2 
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(2) 



-2u„ 



= m_+ 



3u 2 -u, 



1 (1) 



(0) 




00.,= w 2 -3u+u 2 



2 (0) 



-u r u 2 



-3u+u 2 



FIG. 3: Integration ranges for a c i and cr C 2 when U2 > tti (i.e. \u\\ > 1*2])- For a c \ the numbers indicate the fraction value 
in Eq. (1821 and the full lines are where the oji integral is dominant. For a C 2 the numbers in parenthesis indicate the fraction 
value in Eq. 1861 and the arrowed lines are where the u>2 integral is dominant. 



The saddle point of this integral is below — 2u2, hence it is dominated by UJ2 = — 2?i2, i.e. is the same as if the 
latter is also dominated by uo-2 = — 2it2, or less than <r^ if the latter has a saddle point within the integration range. 
Hence cr^ determines the result with 

a cl „ e ~(m+u 2 f/(2A 1+ 2A 2 ) (M2 _ Ui)A2 K _ 2u2 A 1 

a cl ~ e -(«i-«2) 2 /(2A 1 )-2^M 2 (U2 _ Ui)M > _ 2u2 A 1 (85) 
Consider next a c2 which for u>%,u>2 > i s dominated by 

J2 



, U VZJ I VUJ2 v A I C ~ F ZiO . /o/^ 

^ = Z ^ 1 + e t 'l+ tl 2+^l+"2 + e 4ti 1 +2w 1 + g4« 2 + 2^2 /"1^2>° • ' 86 ^ 



The fraction above has values 0, 1, 2 is indicated in Fig. 3 with boundaries shown by the arrowed lines; at these lines 
<9/<9cl>2 is maximal and the corresponding L02 dominate the integral. Hence for lo\ < U2 — u\ the integral is dominated 
by 0J2 = — 2it2 leading to a contribution 
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dcu ie -^^ Al e~ 2u ^ A2 ~ e r 2 <' M . (87) 

The next range is — u\ < u)\ < —u^ — U\ where uji — —uj\ — u\ — 112 dominates, contribution 

U2—U1 

(2) „ I dwie -(^+u 1 +„2) 2 /2A2 e -^/2A 1 _ (88) 



a 



'U2—U1 

This has a maximum within integration range if {v.% — U\)A>2 < —2u 2 A\ with the result 



(2) 



-(« 1+ u 2 ) 2 /(2A 1+ 2A 2 ) (U2 _ Ui)M < ^ 2u2 A 1 (89) 



while if (u2 — Ui)A > —2u 2 A\ the integral is dominated by its lower limit u 2 — u± which is then always smaller then 

(2) 

a c2 ■ In the range —u\ — u 2 < uj\ < — 3ui + U2 the line of maximum L02 = — 3«2 + Ui is at large values of to 2 (see 
Fig. 3) so should give a small contribution; in fact integrating this line even from U2 — u\ yields 

CT (3) ^ f ^ 1 g-("i+"i-3«2) 2 /2A2 e -w 2 /2 J 4 1 ^ e -2« 2 /A 2 -(u2-«l) 2 /2A 1 j-qq) 

J U 2 — Ml 
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where the integrand has a maximum below the integration range and is therefore dominated by the lower integration 
limit. The result in Eq. (|90[l equals also to that of the integrand in (|88[) at its lower limit, hence < af 2 . Finally, 
the range — 3i*i + u 2 < u>i has the line of maximum at w 2 = + 3ui — u%, hence 



(4) 



du> 1 e~ < - UJl+3ui ~~ U2)2 / 2A2 e~ UJ ^ 2Al ~ e -( 3 "i-"2) 2 Mi 



(91) 



— 3«i+u 2 



where again the integral is dominated by its lower limit. This result is smaller then Eq. I|89|) (it is smaller then the 
integrand of (|88|l at ui = —u\ — u 2 , hence the latter is bigger if it has a maximum within integration range). 
Collecting all terms we have, 



C C 2 
C C 2 



max[e- 2 ^, e -(^+^r /(2A 1+ 2A 2 ) 
max[e- 2i1 ', e~^ ul - u ^' M ] 



(u 2 - ui)A 2 < ~2u 2 Ai 
(u 2 - ui)A 2 > -2u 2 Ai . 



(92) 



Eqs. i|85!92l) can be written in terms of a and an anisotropy parameter 77 = Ki/Kq > (For 77 < we note that the 
solutions are symmetric under 77, a — > —77, 1/a). For 77 < 



(Tel ~ Wcl] 16(1+a " T 



a(l-t;)^ 

g c2 — max{ [<j cl ] 



(l + c, + ,,- V a)< ! 

[(Tel] 16 < 1 +°>- } 



(93) 



while for 77 > we have 



(Td 



[(Tel 



(l-a + ri + V c,) z -\-4a(l- V )^ 



a C 2 ~max{[cr c i]" 



(3 + 3^-Q + tlQ) 2 

[(Tel] 16 ' } 



(94) 



Some inspection shows that the latter equation has no solutions (except with a = 0, see below) while for Eq. (|93|l we 
have the following solutions (see Fig. 4): (i) The 2nd term of a c2 (~ er^) identifies a c i exponents and leads to a = 1 
and criticality at a cr = |, i.e. the independent layer solution N = 1. The 2nd term of the a C 2 line is the maximal one 
if 77 < 1 — l/\/2- (h) The 1st term of a C 2 identifies exponents as (1 — 7y) 2 /4 = (1 + a + 77 — ?7a) 2 /16(l + a). This term 
dominates in the a C 2 hue if a < 1, hence the solution 



(95) 



exists for 1 — 1/V2 < 77 < 1/3 with a cr = ;j(l — ??) 2 - (hi) Finally a — is possible, i.e. o c i = and an onset of 
just the k = 7r component a C 2- The solution is then of charges correlated between layers, i.e. the N — 2 rod phase. 



a 



0.5 rj 



FIG. 4: Two layer solutions for the exponent in a C 2 ~ [o"ci] Q in terms of the anisotropy 77. In the range 1 — l/y/2 < 77 < 1/3 
two solutions coexist. 
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FIG. 5: The critical disorder a cr for a two layer system. At 77 < 1 — 1/V2 the transition is to a N = 1 phase at a cr = |. For 
77 > 1 — 1/V2 at u cr = ~(1 — t?) 2 the transition is either to an N = 2 rod phase at 77 > 1/3 or, for 1 - l/y/2 < 77 < 1/3 a mixed 
phase <T C 2 ~ [°"ci] a is possible. At 77 > 1/3 the N = 2 rod solution disintegrates into the N = 1 phase at o-'r' 2 ' = (1 + 77) 2 /16. 

Criticality is at a cr — \(\ — r/) 2 , and from both Eqs. I|93I94I) this solution is valid at all 77 provided it precedes the 
solution (i) with a cr < |, hence 77 > 1 — 1/V2- 

The solutions (i) and (iii) reproduce the energy rationale. We have found here an additional solution (ii) with a 
nontrivial new exponent a in a narrow range 1 — 1/V2 < 77 < 1/3. This solution is a continuous interpolation in a 
between the N = 1 solution (a = 1 at 77 < 1 — l/\/2) and the N = 2 rod solution (a — at 77 > 1/3). Both solutions 
(ii) and (iii) have the same a cr , hence they may be degenerate. 

Solution (iii) allows for an additional phase transition corresponding to the onset of a c \, i.e. the N = 2 rods 
decompose into independent N = 1 charges on each layer. When cr c2 ^ 0, 7*2 and j4 2 are finite, hence the divergent 

terms in the exponent of Eq. (|85|) yield a c \ ~ e - " 1 / 2 " 41 ~ v^i'^ ^ 16<T , hence tTcr' 2 '' = (1 + 7/) 2 /16 allows the onset of 
(T c i at 77 > 1/3 (dashed line in Fig. 5). The energy rationale gives a somewhat higher cr^ 2) = 77 2 /[4(V2-l) 2 ] for this 
N = 2 to N = 1 transition. 

Finally we consider the disorder-temperature phase diagram. The high temperature part of the phase boundary is 
determined by low order renormalization group as disorder is well behaved. Thus, in either Coulomb gas formulation 17 
or in sine-Gordon formulation we find the recursion with scale £ 




t/j 



FIG. 6: Phase diagram for the onset of the N = 1, 2 instabilities for anisotropy 7; = 0.35. At low T two distinct transitions are 
possible, the first being to the rod iV = 2 phase. At high T the independent layer N = 1 transition dominates and eliminates 
the N = 2 phase. 
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The N = 1 solution is determined by one nonzero entry, hence 2 — Kq + <jK^ = 0; for TV = 2 the solution corresponds 
to one nonzero entry per two layers, with the relative sign qp, hence 2 — 2K ± 2^ + 2<jKq = 0. For 77 > the 
dominant transition (i.e. the one at lower temperature) has the upper sign, corresponding to 02c with k — ir. At 
(7 = this has a critical temperature lower than that of N = 1 since K = 1/(1 - 77) < 2 for 77 < \. Therefore the 
range of low a is dominated by the usual N = 1 transition. In Fig. 6 we demonstrate the phase diagram with 77 = 0.35 
where the phases N = 1, 2 compete. 



VI. APPLICATION TO SUPERCONDUCTORS 

A. Layered superconductor without disorder 

The standard model for layered superconductors is the Lawrence Doniach model in terms of the superconducting 
phases on each layer and the electromagnetic vector potential. The latter can be integrated out 9 leading to an effective 
Hamiltonian in terms of pancake vortices, i.e. point singularities in each layer, and a nonsingular Josephson phase. We 
consider here the case without Josephson coupling, where the pancake vortices are not coupled to the Josephson phase. 
If n(r, I) is an integer field of ±1, corresponding to the location of pancake vortices then the vortex Hamiltonian is 9 

h v = Iy, £ n ( r ' ~ r '> 1 - ''M^. + Ec £ n ( r ' ( 9? ) 



2 

r^r' W 



with: 



G ^ = ^irkk) (98) 



2 t hM) a (99) 

where A a h is the magnetic penetration length parallel to the layers and G v (q 7 k) = dj^i J d 2 rG v (r,l)e tkdl+l ' i ' r . The 
core energy is estimated aa 25 i 26 E c « (0.04 — 0.2)r where r = <5>Qd/(47r 2 A 2 h ). 

Note that the k = mode is screened, i.e. G v (q,k) is nonsingular at q — 0. All the other modes are unscreened 
and lead to logarithmic interactions. This is because no screening current can go along z (in the absence of Josephson 
coupling) and thus two pancakes in two different layers cannot screen each others. 

In presence of an external field B along z a flux lattice with a unit cell area a 2 = &o/B is formed. The flux lattice 
is composed of pancake vortices, i.e. point singularities, which are displaced from the 73-th line position R p at the Z-th 
layer into R p + u p ; its Fourier transform is 

u(q, k) =da 2 YYl u l p e l ^ +lkdl (100) 
1 v 

Expanding Eq. l|^7|) to second order in yields the elastic Hamiltonian of the form, 

H e i = \J k J[D L (q, fc)K(q, k)\ 2 + D T {q, k)\u T (q, k)\ 2 ] . (101) 

We will be mainly interested in the case of no Josephson coupling, where the following exact expression holds: 

D L (q, k)P^(q) + D T (q, k)P^(q) (102) 
= ^(q a q f 3G v (q,k) + J2((Q + QhiQ + q)pG v (Q + q,k) - Q a Q fj G v (Q,0))) (103) 

provided we add a short distance cutoff in plane, i.e replace G v (q,k) — > G v (q,k)e~ q ^ I 2 . The conventional elastic 
moduli are then identified as: 

D L (q,k) =q 2 c 11 (q,k) + k 2 z 4 i (q,k) (104) 

D T (q, k) = q 2 c m {q, k) + k 2 z c^ 4 (q, k) (105) 
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where k\ = -^sin 2 (^-). For zero Josephson coupling it is foundSI, where leading terms in q 2 are retained, 

c 66 (<7> k ) = 73— r^2 ( 106 ) 



^2 J 

cu( ff , fc) + c*fe, fc) = ^ 1 + A 2 fc(g2 + ^ 2) (107) 



- ^ rrwri + c - (fc) (108) 



Q^O 

and the last form is in the limit d <C a,X a b- We note that with Josephson coupling the results for C66,cn are 
unchanged, while are modified with a stronger effect on cj 4 . 

We consider first the defect transition in the pure system^. This refers to the proliferation of vacancy interstitial 
pairs (VI), thereby destroying the superconducting order parallel to the layers. These defects correspond to additional 
pancake vortices, denoted by sj(r) on top of the ones forming the flux lattice. These defects couple to the lattice via 
the same coupling of Eq. i|97fl. 

H vac = bi(t)G v (R p + v$-t,1-1') (110) 

To 0-th order in uj 1 the defects feel a periodic potential: 

ni°] c = Yl s i(r)Gv(R p -r,l- I') (111) 

r,p,l,l' 

which fixes the defect position in a unit cell, hence fluctuations of s(q, fc) = dj^i J2 r si(r)e lkdl+lcl ' r involve only q in 
the first Brillouin zone (BZ); in the following (and in Eq. HOlll all q integrals are restricted to the first BZ. Note that 
for vacancies the periodic potential has minima on the flux lines, while for interstitials the minima are in the middle 
of the unit cell. Hence, the core energies of vacancies and interstitials differ, but as they come in pairs, E c refers to an 
average of these core energies. For an isolated pancake vorte:s2f 5 E c ss (0.1 — 0.2)t, while in presence of a flux lattice 
with local relaxation leads to2£ E c w 0.04t. 

Expanding to first order, one finds with the above definitions: 

H vac (s, u) = Hi°l(s) + HW( S , u) + 0(su 2 ) (112) 

The total energy is thus: 

H el (u)+H^ c (s,u)+H v (s) = ^J J D T (q,k)\u T (q,k)\ 2 (114) 

Iff 1 2 

'(D L (q, k)\u L (q, k)\ 2 + -js(q, k)G v (q, k)s(-q, -k) + -^s(q, k)G v (q, k)(-iq)u L (-q, -fc)] (115) 



27 fc 7 9 ^-vi.-v^wv-v ^ -j 1 a 2 d 2 

One can either minimize it to find the (purely longitudinal) deformation of the lattice induced by the defect, 

u vac (q, fc) = iqs(q, k)G v (q, k)/a 2 d 2 D L (q, fc) (116) 

and compute 7i e i + Ti. vac + 7~Lv at the minimum or, since it is Gaussian, simply integrate out the displacements 
U l (q, fc) ■ One finds that the screening of the vortices by the longitudinal displacements of the lattice results in an 
effective interaction energy between the defects: 

«;"(«) = /s(q,t)G:"(,,fe)»(-q,-t) (117) 



-^H'-sMbS) 11181 
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One can connect with the notations of the previous sections (j3 — 1/ksT): 



j(k) = -flimq 2 G e J f (q,k) (119) 



The pure defect transition thus occurs when: 



K ef/ = [ 9(k) = 2 <- T def = i- f G e Jf(k) (120) 



where we define: 



G v (k) = limy G v (q,k)] = f 2 (121) 

G e Jf(k) = km^G^fefc)] = G v (k)(l a , d f^ k) ) (122) 

where we recall k 2 = -^sin 2 ^^). 

It is instructive to consider the "unscreened defect transition" temperature, i.e. formation of pancake vortices in 
the absence of an external field. This is denoted as the vortex transition* with the onset temperature at 

T " = IT I G ^ = o(ft v^ /(2A afc )) (123) 

h(y) = ^ dx- * = 1 - - ? I= (124) 

277 J- 1 + 55% 



In particular for d <C A one has: 



* §d -.r/8. (125) 



2(47rA a6 



The actual superconducting transition is at T c with T v < T c < Tf where Tf is the fluxon transition temperature, 
where Josephson decoupling would occur in the absence of pancake defects^. 

To compare the vortex transition with melting wc use a Lindemann type criterion (with only transverse modes). 

c 2 a 2 = (u 2 ) =T m f f 2 ] = / ln(l + (126) 

Using a circular BZ of volume (27r/a) 2 , hence < < Air/a , 

_ 47TC? , , 47TC? 



A ~ DD " - ^ (87rA ah )2 

A = d / ln (! + ~^%) = d / HI + 1fi 4 ,f°n , n re ) (127) 
J fc c| 4 a 2 fcf J fc 167ra 4 A 2 L)T(0,fc) 

where in the last equation we have used the dispersionlcss value of cqq valid for a d. The scales of the vortex and 
melting transitions are the same, their ratio being T^ e j/T m — A/Airc 2 L . Hence the condition that the defect transition 
occurs before melting and can thus be consistently described is that G%^{q, k) -C G v (q, k). 

Let us now study the true transition with screening. One denotes D L ^ T (k) = D L T (0, k) = k 2 c^ T (0, k) respectively. 
Using the above result, one finds in the q — > limit the exact expressions: 

D T (k) = \{^f ]T(G,(Q,fc)-G,(Q,0))Q 2 (128) 
D L (k) = D T (k) + ~^G v {k) (129) 



One thus has: 



G v {k)-G v {k) DL{k) - ^ 1 + X 2 bk 2 1 + e{k) 

e(k) = a 4 d 2 D T (k)/G v (k) . (131) 
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Hence the condition for Tdef <C T m , which justifies our description of the defect transition, is e(k) <§; 1. We note 
also that for a single 2D layer there is no tilt modulus for the FL and e = 0; hence a 2D FL has VI's at any finite 
temperature. 

Let us first consider the regime 2nd <C a < A relevant for layered superconductors. As shown in the Appendix one 
has in this regime: 

^^'w^*-* 1 11321 

where d = max(d, £/2n) and for d < £ only the first term contributes. This yields: 
e(k) w « 2 ^ + K b kl x A + kllQl nd 2 (l + X 2 ab k 2 z ) 

Note that the relative contribution of the second term becomes significant only for k ~ l/d. The condition that 
e(k) <C 1 for all k is thus met for A/a sufficiently large (high enough field) as: 



X?_ > J_ ln( _?£ 

a 2 An dQo 



->-M — ) (134) 



where c is a constant of order O(l) (which can be estimated from above, with c — 1 when d < £/27r). As long as 
e(fc) <C 1 we find that in all regimes one can estimate: 



* 2 d 2 k 2 & 2 n 2 rl 2 
Gl ff {k) « ^_^_ e(fc) „ J^l n( i + fc 2 /Q 2 ) , 



This yields the estimate of the defect transition for 2nd <C a < A, using Eq. (|f 20(1 at 2nd <C 



(136, 

We use this To as a convenient scale below. We note that (Ilrt5j) is weakly k dependent, hence small anisotropy r/ (Fig. 
1) and the one-layer N=l transition dominates. 

We now estimate the defect transition temperature in the other limit 2nd > a relevant for multilayers. As shown 
in the Appendix one has then: 

v 7 2vra 4 z (l + 2X Q« ) 2 v ; 

This yields: 

<*) = Qode ' Q ° d ^Sly (138) 
The condition for e(k) <C 1 for all k is satisfied when 27rd 3> a. Thus in this limit one finds: 

G eff (k) M #*! fc 2 (139) 



yielding: 



0§ dQ e- d Q° 

Tde/ ~«(TTf^ (140) 



using /, k 2 — 2/d 3, . If in addition A > \/ad one finds: 

def 64tt3A 4 ' 1 j 

The melting criteria Eq. 11271) has now cj 4 which is exponentially small ~ e~®° d , however it enters the logarithm in 
Eq. 1(1271) . i.e. ^4 w 27rd/a is large and T m w c|0oa/(327r 2 A 2 h ), hence T de/ < T m for all 27rd > a. We note that Eq. 
((139(1 implies a significant interlayer coupling with rj (see Fig. 1) close to 0.5; hence disorder favors rod phases at low 
temperatures. 
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B. model with disorder 



We have seen in section IV that disorder can affect a Coulomb gas transition if its correlation diverges at least 
logarithmically with distance. Therefore, disorder that couples directly to pancake defects has a finite correlation and 
has no effect on the defect transition. In particular the vortex transition Eq. (|123[) in the absence of an external 
field is disorder independent. The presence of a flux lattice deformed by point disorder leads to a significant change 
in the disorder as seen by pancake defects. Since each pancake composing the flux lattice is a charge interacting 
logarithmically with the pancake defects, a displaced pancake is equivalent to an addition of +, — charges, i.e. a 
dipole. Hence a disorder deformed flux lattice leads to a quenched dipole disorder seen by the defects, leading to 
logarithmically correlated disorder. 

We consider first disorder within the finite Larkin scale where one can expand in displacement, resulting in a random 
force fi(r) 

H dis = ~J2 [ d 2 rf;(rKG,r) (142) 
i J 

Mr)f l ,(r') = F l _ l ,6 2 (r-T / ) (143) 

where we display only the coupling to the longitudinal component, ul(1,y) (being a suitable continuation of u l p ) and 
/z(r) is the longitudinal disorder component; only the longitudinal mode Ul{1,t) couples to the defects (Eq. I113[) . 
One can write it in Fourier components: 

H dis = ~JJj(k,q)u*(k,q) (144) 



f(k, q)f(k', q') = (2^) 3 <5 2 (q + q')5(k + k')F{k) (145) 

where F(k) = dJ2 t p l6 tkdl '■• no t e that for finite M one has J, = ^jfJ2k an< 4 27T(5(fc + k') = dMSk,k'- It is useful 
to relate F(k) to a previously used*£ dimcnsionlcss disorder parameter s representing point disorder uncorrelated 
between layers. The replicated action has (32ttsTq / da 4 ) J k u a (q, fc)iij*(q, k) with To from Eq. H136(l . Replicating Eq. 

(Q3Ul identifies F(k) = GAnsdTfi /a 4 . 

We now consider the total energy H e i(u) + H vac (s, u) + iJ^ s (u) and determine the u configuration in presence of 
both disorder and defects. The part involving longitudinal displacements reads: 

f 1 1 

Htot = J ^-^pS(q,k)G v (q,k)s(c[,ky+ -^-^s(q,k)G v {q,k)(-iq)u* L (q,k) + 



£l^l\u L (q, k)\ 2 - ~/(q, fcK(q, k) 



and we neglect the random potential seen by the defect itself (which is short range). The relaxed phonon field at the 
minimum energy is: 

u L (q, k) = *q S (q, k) f/^ k) + ^^/(q. *) (146) 

Computing the energy for the defects at the minimum (or equivalently integrating out the displacements) yields the 
same screened interaction Gf,^(q,k) between defects as before and in addition yields the coupling of the vacancy 
to disorder (through the lattice) as in the starting model which allows to identify the correlator A(fc) introduced in 
Section II: 



H vd , ls = -Y,Vi(r) Sl (r) = ~ f V(q,k)s*(q,k) (147) 



An 



V(q, k)V(q', k') - {2ir) 2 6(q + q')(2%)5(k + V)-fA(q, k) (149) 
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Thus in the limit q — > one obtains in general: 



a « - r^^k <151) 

with A(fc) = A(0, k). In the almost fully screened case of interest (2ird <C a < X or 27re? ^> a) we have 
G v {k)/a 4 d 2 D L (k) » 1, hence: 



A(fc) = — (152) 

For usual layered superconductors 27rd <C a we have from Eqs. (|135I136|I for almost all k (k > 1/a) g(fc) = 2(3Tq, 
hence <r e // of Eq. Hfi3[) with iV = 1 becomes <r e // = 4s. Note that cr e // ~ -B, hence defect formation is induced at a 
fixed disorder by increasing the field B. 

We proceed now to study the full disorder problem on all scales allowing for Bragg glass (BrG) properties^. The 
basic assumption is that the long range extra displacement induced in the BrG configuration by the defect is very 
small and one can expand in it. Consider then Hbg{ u ) as the BrG hamiltonian for the u field in presence of disorder 
but in the absence of point defects. We add to it (here u = ul): 

H{u)=H BG {u)+ I h{q,k)u(q,k). (153) 

J g } k 

In particular for the flux lattice problem we identify from Eq. (|113|l 

h fa fc ) = JV^' k ) G v(l, k) ■ (154) 
ercr 

The next order in the displacement expansion is 0(su 2 ) and after integrating out Ui(q, k) leads to s 3 and higher 
order terms; these are neglected in our low density treatment of defects, i.e. large [3E C . 
Then one has the exact (although formal) expansion for the free energy F = —ThiZ: 

F = F BG + [ h(q, k) < «(q, k) > ~ f h(q, k)h(q' , k') < u(q, fe)«(q', k') > c +0{h 3 ) (155) 

Jq,k Z1 Jq,q',k,k' 

where < ... > is thermal average in a particular disorder configuration with no defects and Fbg is the free energy of 
the BG in that configuration and c denotes connected averages; disorder average will follow below. 

In the absence of disorder the second term in Eq. (|155fl is zero and the third one yields the energy which screens 
the initial defect-defect interaction: 



1 f q 2 G v {q,k) 2 2 

|s(q,A:)| ( 156 ) 



iL ' ee " 2d 2 J q ^ k d 2 a i D{q 1 ky 

using < u(q, fc)it(q', k') > c = T/D(q, k), i.e. the screening term in Eq. I|118|l . 

In presence of disorder, the disorder average of the third term in Eq. (|155fl still yields exactly the same screening 
part of the interaction between defects. This is guaranteed by the so-called statistical tilt symmetry of the Bragg 
glass model in the absence of defects, i.e. the statistical invariance of the disorder term in the Hamiltonian under 
u(r, I) — * u(r, I) + 4>(r, I) where <j>(r, I) is an arbitrary function (see e.gi) so that < u(q, k)u{q[, k') > c ~ 5 2 F/6(t> 2 \ < f,=o 
is independent of disorder. 

Since this is an expansion in defect density s(q, k) we can now identify the random potential coupling linearly to 
the defect via the second term (a response of the third term in Eq. (|155fl to defects results in higher order 0(s 3 ) 
terms): 



V(q, k) = --^2 G v{q, k)iq- < u(q, k) > B g (157) 



The correlations are thus (overbar is disorder average): 



47T 



V(q, *)V(q', k>) = {2n) 2 5{q + q'){2%)5(k + k')^-A(q, k) (158) 
A(q, k) = -^±^q*G v (q, k) 2 C BG (q, k) (159) 



23 



where Cbg(<1, k) denotes the disconnected average: 



< u(q, k) >< u(q', k') > = (27r) 2 (Hq + q')(27r)5(fc + k')C BG (q, k) (160) 

At all temperatures except near melting one has:< u >< u > w < uu > as thermal fluctuations are subdominant. 
Therefore we replace the left hand side of Eq. l|ltj(J|) by 1 

C BG (q,k)^l/(q i + q 3 /R c ) (161) 

where q 2 = Cge? 2 + c% 4 k 2 and R c is a Larkin length along c. For q = and large fc > 1/R C , i-e. on short distances 
compared with i? c , this reduces to the previous result Eq. 1)152(1 . while at longer scales the BrG induces interlayer 
disorder correlation as seen by the defects. Replacing l/D 2 L (k) in Eq. (|151|l by Cbg(<1, fc) at g = we obtain (using 
G v (k)/a 4 d 2 D L (k) w 1 as in Eq. 

A « " <162) 

It is instructive to present another derivation of A(fc), valid at T = 0. In general, the disorder potential V(r,l) 
couples to the flux density p(r, u(r, I)) and leads to a Bragg glass configuration UBc(r, I). The addition of a vacancy at 
position R on layer I leads to an energy of U(R) = ^ v J d 2 ru vac (r — R, I' — I)- Vp(r, u BG (r, Z'))F(r, I'). One can now 
see that the force Vp(r, UBg(r, l))V(r, I) has short range correlations. Indeed, at T = we can minimize the disorder 
energy J2i J d 2 rp(r, u(r, l))V(r, I) with the elastic energy Eq. (|ll)lfl to yield uacfr, Z)), hence Vp(r, ubg( t , l))V( r , ~ 
V 2 us G (r, Z)), the latter quantity having clearly short range correlations ~ q 4 /[q 4 + Rj 1 ^ 3 ]- The potential U(R) is 
thus the convolution of a short range correlated random force with the displacement u vac which has a long range 
form: for a single vacancy |u mc (q, k)\ 2 ~ 1/q 2 from Eq. i|lltj|) . Thus one finds that U(R) is logarithmically correlated 
with A(fc) of Eq. i|162[l . Hence the BrG induces an effective disorder correlation between layers on scales longer than 
R c . For weak disorder R c ^> d and the effect in J, A(fc) is negligible, hence the results of the Larkin regime are valid. 

The application of our results to flux lattices depends on the interlayer form of Eq. I|98|) which for a 3> d has the 
form 9 G v \r,l) ~ e~ u l a lnr, i.e. a range of Iq w a/d. For usual layered superconductors 2 with a/d w 10 — 100, we find 
that the AT = 1 phase dominates and cr cr = 1/8. The phase diagram has then the form of Fig. 2 with the magnetic 
field B in the vertical axis 

To achieve N ^ 1 phases the nearest layer coupling should increase. We note that g(k = 0) = since for a straight 
pancake rod the logarithmic interaction is fully screened. Hence J; = and when the range Iq is reduced Jo, Ji 
dominate the sum, i.e. rji — > | when d > a, as in Eq. I|139fl . Direct evaluation of 771 shows that it crosses the 
critical value 1 — l/v2 when d/a w 1, depending weakly on the ratio a/X a b- We therefore propose that flux lattices 
in multilayer superconductors, where d > a can be achieved, may show a rich phase diagram with N > 1 phases. 



VII. DISCUSSION 



We have developed here a variational method and a Cayley tree rationale and applied these to the layered Coulomb 
gas. The variational method is shown to reproduce the defect transition of the single layer as well as demonstrate 
a first order transition within the ordered phase. The latter was so far inferred in the Caylee tree problemSS or in 
the dynamic problemii. To observe this transition one needs to induce defects in the system, e.g. by finite size or 
dynamics. We also show that this line survives in the disordered phase, showing a crossover in the defect density 
dependence on temperature or disorder. 

The variational scheme has been extended to two layers, confirming essentially the energy rationale. Near the 
onset to the N = 2 rod phase we find in a narrow interval a curious phase with a new exponent relating the two 
components of the order parameter. We consider then the variational scheme as reliable for the main features of the 
phase diagram, i.e. the sequence of transitions into rod phases (Fig. 1). 

Our results are relevant to flux lattices where we find the phase boundaries and propose that for 2nd > a new N > 1 
phases can be manifested. Our derivation assumes (i) dislocations are neglected, and (ii) the Josephson coupling is 
neglected. Assumption (i) implies that the melting transition is at higher temperature or disorder than those of 
the defect transition. This has been justified for the pure case in section VIA showing that Td e f <C T m if either 
2nd <C a < X or 2nd ^> a. We assume that the same holds for disorder induced melting, though the latter is less 
understood. 

We discuss next assumption (ii), i.e. the effect of the interlayer Josephson coupling J. In the absence of VI a layer 
decoupling was founc&£ where J vanishes on long scales. At this transition the width of a Josephson flux line diverges 
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and its fluctuations renormalize J to zero. A complete description should allow for both VI defects and Josephson 
vortex loops which would combine to form 3-dimcnsional defect loops. We expect then that the defect and decoupling 
transitions merge into one transition at T c , above which both the renormalized J is zero as well as a finite VI density 
rid appears. 

In fact a transition to a " supersolid" phase in a flux lattice in isotropic superconductors was proposed^ where a 
finite density of defect loops proliferate and a related " quartet" dislocation scenario was suggested^. In the supersolid 
desription a finite line energy competes with the entropy of the wandering line, both being linear in the defect length. 
The resulting transition temperature is comparable to that of meltingifi, hence it is uncertain if this scenario is 
possible. 

In our VI transition the competing energies and entropies are logarithmic in the VI separation, rather than linear. 
If a Josephson coupling is added, naively a linear term is added since a flux line connecting the VI pair is formed. 
However, near decoupling the renormalized J varies as a power of scale, hence we expect that the free energy of a 
flux loop to be nonlinear in size, modifying significantly the supersolid transition at least in the small J case. We also 
show now that, in contrast with the supersolid scenario, T c can be well below melting. 

We note first that in the pure system the decoupling transition is at Tdec — 8Td e f (for d <C a <C A) while its critical 
disorder (at T = 0) is al£ a dec = 2 = Wag e f, hence the a — T boundary of the defect transition is below that of 
decoupling in both the a, T coordinates. The disorder-temperature "phase diagram" has therefore 3 regions, separated 
by the two lines Td e fis T ) and Tdecis 7 )'- (i) decoupled and defected phase at high T or high a, (ii) between the lines 
Tdef {.<?)■> Tdec{ cr )> an d (hi) a coupled defectless phase at small T and small a. This "phase diagram" is inconsistent in 
the sense that Tdef if) is derived in the absence of J, while Td ec (cr) is derived in the absence of VI defects. 

We show next that T^/ < T c < Td ec - In phase (i) J — > and rid is relevant in the RG sense. This is a consistent 
description since J = is assumed in the VI description, hence region (i) is a disordered phase. In region (hi) rid ~ *■ 
while J is relevant, again a consistent scenario since J being relevant is shown assuming rid — 0. However, in region 
(ii) both rid and J are relevant, hence seperate "decoupling" and "defect" descriptions are inconsistent and a single 
combined transition within region (ii) is expected, i.e Tdef < T c < Tdec- Since both Tdef, Tdec are well below melting 
for a <C A, we conclude that T c is also well below melting. 

In fact we can estimate T c by an argument as used in the B — case 9 . Consider the VI correlation length 
(d Ri n, for J — (which diverges at Tdef) and the Josephson correlation length £j (which diverges at Tdec)- 
Consider a temperature for which £j < £<j; £j is the scale at which J/T is renormalized to strong coupling w 1, e.g. 
in 1st order RG— 

O « a(T/J) 1 ^ 1 - T / T ^ . (163) 

The Josephson term J cos(8 ns + 8 S ) involves both the nonsingular phase ns and the singular one S due to VI pairs. 
If £j < £d VI pairs are not seen on the scale between a and £j, renormalization of J cos(9 ns ) can proceed till strong 
coupling is achieved, i.e. the phase is ordered. If instead £j > £j VI defects interfere in the J renormalization and 
disorder the system. Hence T c is estimated by £j w From Eq. I)55ll for the pure case, 

£ d W a ( e /^=)l/[2(l-WT)l ( ( 164 ) 

hence T c is near Tdef if J is sufficiently small, 

J < Te- E " /T . (165) 

For Bi 2 Sr 2 CaCu 2 O s we estimate^ J w 0.1K, E c w 10 3 K which for relevant T = 10 - 100K does not satisfy 
Eq. ^ 165(1 . i.e. the transition is near Tdec- However, for multilayers such as (Bi2Sr2CaCu 2 0%) m (Bi 2 Sr 2 CuO^) n ^ 
the semiconducting layers of Bi 2 Sr 2 CuO§ reduce J by a factor e~~ d ^ where d ~ n includes now the thickness of 
the semiconducting layers. Hence, for a few such layers the condition (|165[) is already satisfied and T c is near Tdef- 
Therefore, multilayer systems are excellent candidates for observing the VI transition with interlayer defects being 
either uncorrelated, when l-nd < a, or in correlated N > 1 rod phases, when 2ird > a. The latter condition is in fact 
easier to realize in these multilayers where d is larger. Increasing d too much will push down the coupled phase to 
very low temperatures, hence the optimal case for study are multilayers with 2nd ~ a. 

Acknowledgements: This research was supported by THE ISRAEL SCIENCE FOUNDATION founded by the Israel 
Academy of Sciences and Humanities. B. H. thanks the Ecole Normale Superieure for hospitality and for support as 
visiting professor during part of this work. 



25 



APPENDIX A: EVALUATION OF SOME QUANTITIES 

Let us estimate: 

D T (k) = \{^f ( G -W' fc ) - G -(Q> °))Q 2 ( A1 ) 

in various regimes. Implicit in the sum is a cutoff at large Q w 27r/£. One has: 

^^'^S^TTTky-TTTb) 1 (A2 » 



= ^ V f - - (A3) 

8?rA 2 a 4 ^ 1 J_ d sinh(Qrf) 1 rf sinh(Qrf) x 

Q^O 1 H^Q sinh^Od^+sin^/M^ 1 ) 1 I" H^Q sinh 2 fOd/2) / 



Let us first consider the case d < £o> d/a <C 1. Then it simplifies into: 

-A 2 Q 2 A 2 (^ 
+ A 2 Q 2 + l + A 2 (Q 2 + fc 2 ) 



UT{k) ~ SttA 2 ^ 2-, 4 + A 2 Q 2 + 1 + A 2 (Q 2 + A: 2 ) j (A4) 



If we now further consider the case where A > a it becomes: 

a 2 . — i 1 n . r a /t 



Q#0 

Thus wc find, for (i<(o,(f/«<l and A > a that: 



DT{k) M 32^W ln( l + (6fc z /2.) 2 ) (A6) 
where Qo = 2-7r/a is the lowest term in the Q sum. In general for d > £ one has: 

n (/:) ~ *§ ^ -^Q 2 A 2 (Q 2 + fc 2 ) ^ 1 i 

U 8,A 2 a4 ( 4 4 + A 2 Q 2 + 1 + A 2 (Q 2 + m )+ 2^ 1 + sfo 1+4e ^L^ /2) l + 53 fc ( j 

The first term can be estimated for A > a and the second with no assumption: 

Using the previous calculation this yields: 

and one can compute $(fc) in two limits: 
In the case 2nd <C a one finds: 

^'-sg^^l^l (A10) 

r x u 2 

F[ X ,y] = j i ^(T+^e- (AH) 
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Since A ^> d seems natural in that case one gets: 



$(fc) 



8Q, 



(A12) 



In the opposite case 2nd a the sum is dominated by the two shortest Q of length Q 0l hence 



D T {k) « $(fc) 




1 



(A13) 
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